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[i]  The  dynamical  core  of  the  Naval  Research  Laboratory  (NRL)  Spectral  Element 
Atmospheric  Model  (NSEAM)  is  coupled  with  full  physics  and  used  to  investigate  the 
organization  and  propagation  of  equatorial  atmospheric  waves  under  the  aquaplanet 
conditions.  The  sensitivity  of  the  model  simulation  to  the  amount  of  horizontal  viscosity, 
distribution  of  the  vertical  levels,  and  selected  details  of  the  precipitation  physics  is 
examined  and  discussed  mainly  utilizing  simulated  convective  precipitation  with  the  aid  of 
time-longitude  plots  and  the  spectral  diagrams  designed  by  Wheeler  and  Kiladis  ( 1 999).  It  is 
shown  that  the  simulation  of  the  Kelvin  wave  and  Madden-Julian  Oscillation  depends 
strongly  on  the  details  of  the  vertical  level  distribution  and  the  choice  of  parameters  in  the 
convective  parameterization.  Efforts  are  made  to  calibrate  the  new  model  to  capture  the 
essential  interaction  between  the  dynamics  and  physics  of  the  atmosphere.  The  speed  and 
spectrum  of  the  eastward  propagating  Kelvin  waves  and  the  signature  of  the  Madden-Julian 
Oscillation  simulated  by  the  new  model  reveal  main  features  similar  to  those  predicted  by 
the  simplified  theory  and  found  in  limited  observations.  This  study  attempts  to  understand 
the  significant  variability  found  among  the  aquaplanet  simulations  by  various  global 
atmospheric  models  and  highlights  the  uncertainties  concerning  convective  processes  and 
their  coupling  to  large-scale  wave  motion  in  large-scale  models  of  the  atmosphere. 
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1.  Introduction 

[2]  The  U.S.  Navy  has  developed  a  new  high-accuracy 
global  atmospheric  model,  which  scales  efficiently  on 
current  and  future  state-of-the-art  computing  platforms: 
the  Naval  Research  Laboratory  (NRL)  Spectral  Element 
Atmospheric  Model,  or  NSEAM.  Its  dynamical  core 
[Giraldo  and  Rosmond,  2004;  Giraldo,  2005]  is  based 
on  the  “spectral  element”  method  projected  to  three- 
dimensional  Cartesian  coordinates,  which  effectively  elim¬ 
inates  the  pole  singularity  problem  of  spherical  coordinates 
and  combines  the  local  domain  decomposition  property  of 
finite  element  methods  with  the  high-order  accuracy  of 
spectral  transform  methods. 

[3]  An  important  advantage  of  spectral  element  (SE) 
models  is  that  the  solution  of  the  global  matrix  problem 
required  by  the  semi-implicit  method  is  straightforward  and 
computationally  efficient  for  either  the  hydrostatic  or  non¬ 
hydrostatic  equations.  The  semi-implicit  method  for  either 
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an  SE  or  spherical  harmonics  (SH)  model  requires  the 
solution  of  a  three-dimensional  Helmholtz  operator.  How¬ 
ever,  the  only  reason  why  SH  models  are  competitive  with 
SE  models  (or  grid  point  models  for  that  matter)  for  the 
hydrostatic  equations  is  that  SH  models  do  not  require  the 
solution  of  a  matrix  problem  (the  Helmholtz  operator  is 
solved  exactly  for  constant  coefficients).  Still,  SE/grid  point 
models  are  competitive  with  SH  models  at  high  resolution 
and  on  distributed  memory  architectures.  In  the  nonhydro¬ 
static  case,  the  SH  models  now  require  the  inversion  of  a 
matrix  which  adds  additional  computational  cost  to  the 
model.  The  SE  models,  on  the  other  hand,  already  include 
this  overhead  and  thus  transitioning  to  a  nonhydro  static 
system  incurs  only  minimal  additional  computational  cost. 

[4]  NSEAM  can  adopt  any  horizontal  model  grid,  fixed 
or  variable,  and  various  time  integrators  such  as  Eulerian 
[Giraldo  and  Rosmond,  2004],  semi-implicit  [Giraldo, 
2005],  or  hybrid  Eulerian-Lagrangian  semi-implicit 
[Giraldo,  2006]  methods.  Its  spectral  element  formulation 
maintains  the  high-order  accuracy  of  spherical  harmonics 
that  the  current  Navy’s  global  atmospheric  model  is  based 
on,  while  it  offers  flexibility  to  employ  any  form  of  variable 
grid  to  enhance  horizontal  grid  resolutions  in  strategic 
regions.  Its  dynamical  core  scales  efficiently,  i.e.,  allows 
the  use  of  large  numbers  of  processors,  and  was  validated 
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Figure  1.  An  orthographic  view  of  the  hexahedral  grid 
points  used  in  this  study  for  NSEAM.  Six  elements  on  one 
face  in  one  direction  and  8th  polynomial  order  of  the  basis 
functions  are  selected  for  this  study.  This  horizontal 
resolution  is  roughly  equivalent  to  the  triangular  spectral 
truncation  at  wave  number  54  (i.e.,  T54)  or  the  grid  point 
resolution  of  2.2°. 


using  various  barotropic  and  baroclinic  test  cases  [Giraldo 
and  Rosmond,  2004;  Giraldo,  2005],  The  model  can  be 
discretized  vertically  with  any  grid,  but  the  mass  and  energy 
conserving  flux  form  of  the  finite  difference  method  on  the 
terrain-following  (a)  coordinate  is  first  selected  for  compar¬ 
ison  with  existing  models. 

[5]  For  this  study,  we  modify  and  expand  the  NSEAM 
dynamical  core  to  include  the  physics  package  utilized  in 
the  operational  version  of  the  Navy  Operational  Global 
Atmospheric  Prediction  System  or  NOGAPS  [ Hogan  and 
Rosmond,  1991],  but  without  the  land  surface  parameteri¬ 
zation  (for  recent  model  physics  improvements  refer  to 
Peng  et  al.  [2004],  Hogan  and  Pauley  [2007]  and  Kim 
[2007]).  We  select  the  hexahedral  grid  that  consists  of  six 
faces  of  a  hexahedron,  each  of  which  contains  a  desired 
number  of  quadrilateral  elements  (see  Figure  1)  following 
Giraldo  and  Rosmond  [2004]. 

[6]  Among  various  methods  to  validate  an  atmospheric 
model  is  to  force  the  model  under  controlled  and  simplified 
sets  of  boundary  and  initial  conditions  so  that  the  results  can 
be  interpreted  in  a  relatively  straightforward  manner  and 
also  intercomparable  to  those  of  other  similar  models, 
although  the  correct  solution  is  not  quite  known  except  by 
simplified  theory  and  limited  indirect  observations.  A  good 
recent  example  is  the  aquaplanet  experiments  proposed  by 
Neale  and  Hoskins  [2001a,  2001b]  in  which  the  Earth  is 
covered  with  water  only.  These  experiments  provide  a 
useful  and  convenient  test  bed  for  investigating  the  inter¬ 
action  between  the  dynamics  and  physics  in  atmospheric 
models. 

[7]  In  this  study,  we  validate  NSEAM  by  configuring  it 
for  aquaplanet  experiments  mainly  following  Neale  and 


Hoskins  [2001a],  We  perform  various  sensitivity  experi¬ 
ments  in  order  to  understand  and  improve  the  aquaplanet 
simulation.  Sensitivity  of  aquaplanet  simulations  to  hori¬ 
zontal  resolution  was  studied  previously.  For  example, 
Lorant  and  Royer  [2001]  found  from  general  circulation 
model  (GCM)  experiments  that  with  higher  resolution  the 
convective  cells  are  more  intensified  and  concentrated, 
being  accompanied  by  improved  simulation  of  equatorial 
waves  that  modulate  near-equatorial  convection.  Sensitivity 
of  aquaplanet  simulations  to  vertical  resolution  was  also 
studied  earlier.  For  instance,  Inness  et  al.  [2001]  compared 
between  19  and  30  (unevenly  spaced)  layer  versions  of  their 
GCM  and  discussed  its  implications  for  Madden-Julian 
Oscillation  (MJO)  in  view  of  the  moisture  budget.  They 
reported  that  the  effect  of  convection  is  to  moisten/dry  the 
lower  troposphere  in  their  19/30  layer  simulations,  respec¬ 
tively.  In  the  present  study,  we  investigate  further  sensitivity 
of  the  aquaplanet  simulations  to  selected  details  of  the 
model  such  as  the  vertical  distribution,  as  well  as  the 
resolution,  of  the  model  levels  and  the  way  the  lifting 
condensation  level  (LCL)  is  calculated  or  utilized  in  the 
context  of  the  frequency  and  propagation  of  simulated 
Kelvin  waves  and  MJO. 

[8]  In  section  2,  we  outline  the  coupling  of  the  NSEAM 
dynamical  core  with  the  NOGAPS  physics  including  the 
new  horizontal  viscosity  that  replaces  the  original  horizontal 
diffusion.  In  section  3,  we  describe  the  setup  of  the  aqua¬ 
planet  experiments  and  present  the  results  of  the  sensitivity 
experiments  in  terms  of  the  viscosity,  vertical  distribution  of 
the  model  levels  and  some  details  of  the  precipitation 
physics  relevant  to  LCL.  We  also  discuss  our  efforts  to 
simulate  the  Kelvin  waves  and  MJO  in  this  section.  We  end 
the  discussion  in  section  4  by  giving  concluding  remarks 
and  a  short  summary  of  additional  sensitivity  experiments 
perfonned  but  not  presented  in  this  study.  Appendix  A 
includes  the  derivation  of  the  viscosity  operator. 

2.  NSEAM  With  Physics 

[9]  The  coupling  of  physics  with  NSEAM  dynamics  is 
handled  by  a  fractional  step  approach,  as  similarly  done  in 
NOGAPS.  The  process  in  continuous  space  can  be  denoted  as 

%  =  SD{q)+SP{q)  (1) 

where  q  is  any  prognostic  physical  variable,  S  represents  the 
forcing,  and  D  and  P  denote  dynamics  and  physics, 
respectively.  We  first  solve 

(2) 

and  then  use  this  new,  auxiliary  solution  to  solve 

|  =  W)  (3) 

which  represents  the  forcing  by  physical  parameterizations 

and  is  evaluated  in  the  fine  mesh,  i.e.,  at  the  Legendre- 

Gauss-Lobatto  (LGL)  points. 
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Figure  2.  The  vertical  grids  used  for  the  NSEAM 
aquaplanet  simulations,  (a)  20  evenly  spaced  a  levels,  (b) 
30  evenly  spaced  a  levels,  and  (c)  30  unevenly  spaced  a 
levels  with  tightly  spaced  levels  near  the  surface  and  top  and 
loosely  spaced  intermediate  levels. 


[10]  We  formulate  a  high-order  hyperviscosity  for  physics 
by  generalizing  the  low-order  approach  introduced  by 
Giraldo  [1999]  as 

^[=S(q)+^2Kq  (4) 

where  K  is  the  order  of  the  hyperviscosity  operator  (in  this 
study  we  choose  K  =  1,  i.e.,  second-order),  and  n  =  1  x  106 
(m2  s-1)  is  the  default  viscosity  value  used  in  this  study 
(the  choice  of  this  value  will  be  explained  in  section  3.2). 
The  expressions  for  the  viscosity  operator  in  equation  (4) 
are  derived  in  Appendix  A. 

[n]  For  the  horizontal  grid,  we  choose  6  elements  on  one 
face  in  one  direction  and  8th-order  polynomial  basis  func¬ 
tions  (Figure  1).  The  average  horizontal  resolution  of  the 
grid  is  roughly  equivalent  to  the  triangular  spectral  trunca¬ 
tion  at  wave  number  54  (i.e.,  T54)  or  the  grid  point 
resolution  of  2.2°.  We  use  20  or  30  vertical  levels,  either 
evenly  spaced  or  unevenly  spaced  in  a  with  the  model  top 
pressure  of  1  hPa  (Figure  2).  The  time  integration  is  done 
semi-implicitly  with  a  second-order  backward  difference 
method.  The  default  time  step  is  300  s  while  the  radiation 
parameterization  is  called  every  hour. 

3.  Aquaplanet  Experiments 

3.1.  Experimental  Setup 

[12]  NSEAM  is  configured  to  perform  aquaplanet  simu¬ 
lations  basically  following  the  Aquaplanet  Experiment 
(APE)  Intercomparison  Project  recommendations  [Neale 
and  Hoskins,  2001a,  2001b]:  The  Earth  is  covered  with 
water  with  no  orography,  land  or  sea  ice.  The  sea  surface 
temperature  prescribed  is  zonally  uniform  and  symmetric 
with  respect  to  the  equator  (Figure  3)  as  in  the  “control” 
case  of  Neale  and  Hoskins  [2001a],  Earth  eccentricity  and 
obliquity  are  set  to  zero.  The  solar  insolation  is  fixed  to  the 
equinoctial  condition  with  the  solar  constant  of  1365  WnT2. 
We  include  the  diurnal  cycle  in  this  study.  C02  is  prescribed 
to  the  amount  of  348  ppmv  following  the  Atmospheric 
Model  Intercomparison  Project  (AMIP  [Gates,  1992])  II 
[see  also  Weave,  2004].  The  ozone  is  also  taken  from  AMIP 
II  and  zonally  averaged,  but  is  symmetrized  with  respect  to 
the  equator  for  this  study. 

[13]  Furthermore,  we  take  the  year  2005  averages  of  the 
humidity  and  air  temperature  from  the  National  Centers  for 
Environmental  Prediction/National  Center  for  Atmospheric 
Research  (NCEP/NCAR)  Reanalysis  data  [Kalnay  et  al., 
1996],  which  for  this  study  are  also  zonally  averaged  and 
symmetrized  with  respect  to  the  equator.  We  impose  this 
latitudinal  symmetry  of  the  ozone,  humidity  and  air  tem¬ 
perature  to  remove  the  initial  asymmetry  so  that  any 
asymmetry  that  develops  during  simulation  is  due  solely 
to  the  model  physics.  The  simulations  are  run  for  1  year 
with  the  output  saved  every  12  h.  The  first  30  days  are 
considered  a  spin-up  period  and  excluded  in  the  averages  of 
the  12-h  interval  outputs. 

3.2.  Sensitivity  to  Horizontal  Viscosity 

[14]  Horizontal  viscosity  is  a  crude  way  of  parameterizing 
the  horizontal  effects  of  turbulent  mixing  as  well  as  sup¬ 
pressing  computational  noise.  For  idealized  test  cases  with 
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Figure  3.  (a)  The  sea  surface  temperature  (10  1  K)  prescribed  for  the  aquaplanet  experiments  with 
NSEAM,  proposed  for  the  Aquaplanet  Experiment  (APE)  Intercomparison  Project  by  Neale  and  Hoskins 
[2001a].  (b)  The  zonal  mean  profile  of  the  temperature  with  respect  to  the  latitude.  The  maximum  is 
300  K  (27°C)  over  the  equator  and  the  minimum  is  273  K  (0°C)  at  the  poles. 


NSEAM,  Giraldo  [2005]  used  the  viscosity  coefficient  of 
/x  =  7  x  105  m2  s_1,  based  on  the  suggestion  by  Polvani  et 
al.  [2004]  for  baroclinic  instability  waves.  In  this  study,  we 
experiment  a  range  of  viscosity  values  with  its  maximum 
guided  by  this  value.  We  perform  a  series  of  sensitivity 
experiments  to  determine  the  optimal  magnitude  of  the 


viscosity  given  by  equation  (4)  with  the  Laplacian  operator 
expressed  by  equation  (A2)  in  Appendix  A. 

[is]  Figure  4  presents  the  time-longitude  (or  Hovmoeller) 
plots  of  the  convective  precipitation  averaged  over  the 
equator  (5°S-5°N)  for  the  first  90  days  of  the  simulations 
with  twenty  evenly  spaced  vertical  levels.  We  choose  a 
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Figure  4.  Time-longitude  (or  Hovmoeller)  plot  from  days  0  to  90  of  the  convective  precipitation 
(10  1  mm/d)  averaged  for  an  equatorial  band  (5°S-5°N)  obtained  from  the  NSEAM  aquaplanet 
experiments  with  20  evenly  spaced  vertical  levels  using  the  viscosity  value  (/r:  m2  s-1)  of  (a)  0,  (b)  1  x  1 05, 
(c)  5  x  105,  and  (d)  1  x  106. 
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Figure  5.  The  convective  precipitation  (10  1  mm/d)  averaged  over  days  30-365  from  the  NSEAM 
aquaplanet  experiments  with  20  evenly  spaced  levels  using  (a)  no  horizontal  viscosity  and  (b)  second- 
order  viscosity  of  p  =  1  x  106  (m2  s"  ').  Superimposed  are  the  points  of  the  spectral  element  grid. 


90-day  interval  for  comparison  with  Neale  and  Hoskins 
[2001b,  Figure  3a],  It  is  seen  that  after  about  20  days  the 
simulations  become  quite  steady  (not  shown,  but  at  least  for 
2  years).  Without  the  viscosity  (Figure  4a),  the  organization 
and  grouping  of  precipitating  convective  clouds  are  not  so 
pronounced.  A  careful  look  at  the  details,  however,  reveals 
the  individual  clouds  (or  very  small  clusters)  propagating 
westward  (upward  and  to  the  right  in  Figure  4a)  whereas  the 
loosely  grouped  clusters  propagate  eastward  (downward 
and  to  the  right  in  Figure  4a). 

[i6]  When  a  small  amount  of  viscosity  is  added  (Figure  4b; 
p  =  1  x  105  m2  s~‘)  the  clustering  and  propagation  become 
somewhat  more  noticeable.  With  greater  viscosity,  the 
eastward  propagation  of  the  cloud  clusters  is  quite  evident 
while  the  westward  propagation  on  the  contrary  is  becoming 
less  visible  as  compared  in  Figure  4c  (p  =  5  x  105  m2  s"1) 
and  Figure  4d  (p  =  1  x  106  m2  s-1).  With  the  larger  values 


(Figures  4c  and  4d),  the  organization  and  propagation 
features  are  similar  to  those  of  other  studies  [e.g.,  Neale 
and  Hoskins,  2001b],  Throughout  the  remainder  of  this 
study,  we  choose  p  =  1  x  106  m2  s_I  as  the  default  value 
for  the  viscosity  since  it  produces  more  organized  patterns 
in  the  diagnostics  that  will  be  presented  later  in  this  study. 

[17]  Furthermore,  the  SE  grids  suffer  from  a  well-known 
problem  that  the  horizontal  grid  structure  becomes  visible 
by  concentrated  values  of  a  physical  variable,  notably  the 
precipitation,  due  to  the  reduced  grid  intervals  along  the 
element  boundaries.  The  viscosity  can  alleviate  this  “grid 
imprinting”  phenomenon  as  shown  in  Figure  5,  which 
compares  the  convective  precipitation  fields  averaged  over 
days  30  through  365  in  an  equatorial  band  simulated 
without  (Figure  5  a)  and  with  (Figure  5b)  using  the  default 
value  of  the  viscosity.  Similar  results  were  found  from  the 
NCAR’s  Community  Atmospheric  Model  (CAM)  with  an 
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Figure  6.  The  convective  precipitation  ( 1 0  mm/d)  averaged  over  days  30-365  from  the  NSEAM 
aquaplanet  experiments  with  (a)  20  evenly  spaced  (L20e),  (b)  30  evenly  spaced  (L30e),  and  (c) 
30  unevenly  spaced  vertical  levels  (L30u)  given  in  Figure  2. 


SE  dynamical  core  (M.  Taylor,  personal  communication, 
2007).  These  and  other  (not  shown)  sensitivity  experiments 
suggest  that  the  viscosity  plays  an  important  role  in  orga¬ 
nizing  convection. 

3.3.  Sensitivity  to  Vertical  Level  Distribution 

[is]  We  now  investigate  the  effect  of  the  vertical  structure 
of  the  model.  Figure  6  compares  simulated  convective 


precipitations  from  the  experiments  with  the  three  distribu¬ 
tions  of  the  vertical  grid  introduced  in  Figure  2,  which  were 
averaged  for  about  a  year  except  for  the  startup  period  of  the 
first  30  days  (i.e.,  averaged  over  days  30  through  365).  The 
evenly  spaced  20  and  30  vertical  level  cases  (L20e  and 
L30e,  hereafter;  Figures  6a  and  6b,  respectively)  reveal 
more  concentrated  precipitation  over  the  equator  than  the 
unevenly  spaced  30-level  case  (L30u,  hereafter;  Figure  6c), 
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Average  over  700-950  hPa 


Figure  7.  (a)  The  difference  between  L30e  and  L20e  of 
the  water  vapor  (10-1  g/kg)  and  (b)  the  comparison  between 
L30e  and  L20e  of  the  cumulus-induced  heating  rate 
(10  1  K/d)  derived  from  the  temperature  tendency  by  the 
Emanuel  convective  cloud  parameterization.  The  data  have 
been  averaged  between  700  and  950  hPa  over  days  30-365 
of  the  NSEAM  aquaplanet  simulations. 

and  the  maximum  precipitation  is  the  greatest  in  L30e  and 
the  smallest  in  L30u. 

[19]  We  first  look  into  the  difference  between  L20e  and 
L30e  (Figures  6a  and  6b).  With  more  model  levels  (and  thus 
layers)  in  L30e,  there  is  more  precipitation  along  the  narrow 
band  over  the  equator.  Figure  7a  shows  the  difference 
between  L30e  and  L20e  in  the  amount  of  moisture  present 
in  the  lower  tropospheric  model  layers  (averaged  over  700- 
950  hPa  that  includes  the  LCL).  The  moisture  over  the 
equator  in  the  boundary  layer  is  greater  in  L20e  than  in 
L30e.  This  result  is  consistent  with  the  finding  of  lnness  et 
al.  [2001]  in  that  more  vertical  model  levels  (i.e.,  smaller 
vertical  grid  intervals)  induces  relative  drying  of  the  layers 
due  to  convection.  Whereas  the  amount  of  moisture  itself 
over  the  equator  is  less  in  L30e  (Figure  7a),  the  moisture 
convergence  at  850  hPa  (Figure  8;  represented  by  negative 
values  in  bluish  colors)  is  greater  in  L30e  (Figure  8b)  than 
in  L20e  (Figure  8a).  Similar  results  are  found  at  nearby 
vertical  levels  (not  shown).  These  results  imply  stronger 


vertical  moisture  transfer  in  L30e  induced  by  greater  con¬ 
vergence  (Figure  8b)  along  the  equator.  Therefore,  the  L30e 
grid  provides  more  favorable  condition  than  L20e  grid  for 
the  saturation  and  thus  condensation,  leading  to  greater 
precipitation  over  the  equatorial  band  (Figure  6b). 
Figure  7b  compares  the  zonally  averaged  heating  rate 
induced  by  cumulus  (also  averaged  over  700-950  hPa  for 
days  30-365),  which  is  derived  from  the  temperature 
tendency  calculated  by  the  Emanuel  convective  cloud 
parameterization  employed  in  our  model  [ Emanuel  and 
Zivkovic-Rothman,  1999;  Peng  et  al.,  2004],  The  cumu¬ 
lus-induced  heating  is  significantly  greater  in  L30e  over  the 
equatorial  band,  which  is  due  to  the  greater  condensational 
heating  associated  with  the  stronger  precipitation  in  L30e 
than  in  L20e. 

[20]  Figures  9a  and  9b  show  the  Hovmoeller  plots  of  the 
convective  precipitation  from  L20e  and  L30e,  respectively, 
averaged  over  the  equator  (5°S-5°N)  from  simulation  day 
30  to  120  (now  excluding  the  startup  period  of  the  first 
30  days).  There  is  a  clear  difference  between  L20e 
(Figure  9a)  and  L30e  (Figure  9b)  in  the  speed  of  the 
eastward  propagating  waves  (i.e.,  Kelvin  waves)  while  the 
precipitation  is  of  similar  order  of  magnitude.  With  L20e  it 
takes  slightly  more  than  two  months  on  average  for  the 
Kelvin  waves  to  circle  the  Earth  along  the  equator,  but  with 
L30e  it  takes  more  than  four  months.  This  speed  difference 
can  be  explained  by  an  argument  involving  the  vertical  scale 
of  the  convective  systems  as  represented  by  the  equivalent 
depth  and  associated  vertical  modes  under  the  shallow  water 
theory.  We  calculate  the  vertical  percentage  distribution  of 
convective  cloud  top  pressures  (Figure  10)  as  done  by 
lnness  et  al.  [2001,  Figure  15]  except  that  our  calculation 
is  from  5°S  to  5°N  and  for  every  12  h.  The  comparison 
between  L20e  and  L30e  confirms  that  L20e’s  major  cloud 
tops  (at  300  hPa)  are  roughly  twice  deeper  than  L30e’s  (at 
500  hPa),  which  implies  about  twice  deeper  vertical  scale 
and  higher  propagation  speed  in  L20e,  consistently  with 
Figures  9a  and  9b.  The  propagation  speeds  of  Kelvin  waves 
in  both  L20e  and  L30e  are,  however,  fairly  low  com¬ 
pared  with  the  observed  atmospheric  Kelvin  waves  [e.g., 
Masunaga  et  al.,  2006]  and  other  simulations  [e.g., 
Neale  and  Hoskins,  2001b;  Moncrieff  et  al.,  2007]. 

[21]  Next,  we  investigate  the  unevenly  spaced  30-level 
case  (i.e.,  L30u)  in  comparison  with  the  evenly  spaced 
counterpart  (i.e.,  L30e).  The  moisture  over  the  equator  in 
the  boundary  layer  (averaged  over  700-950  hPa  for  days  SO¬ 
SOS)  is  greater  in  L30u  than  in  L30e  (Figure  11a),  meaning 
more  moisture  is  present  in  the  boundary  layer  in  L30u. 
However,  the  precipitation  pattern  of  L30u  case  (Figure  6c) 
is  weaker,  while  broader,  than  L30e  (Figure  6b)  over  the 
equatorial  band.  Rather,  the  difference  in  the  precipitation  is 
consistent  with  the  difference  in  the  moisture  divergence 
between  L30e  (Figure  8b)  and  L30u  (Figure  8c);  that  is, 
much  weaker  and  broader  precipitation  in  L30u  is  attribut¬ 
able  to  the  much  weaker  and  broader  moisture  convergence 
in  the  boundary  layer.  Moreover,  Figure  lib  shows  that  the 
cumulus-induced  heating  is  greater  in  L30e  than  in  L30u 
(also  averaged  over  700-950  hPa  for  days  30-365).  This 
result  is  again  consistent  with  L30e’s  stronger  low-level 
moisture  convergence  (Figure  8b)  and  precipitation 
(Figure  6b)  along  the  equator  in  view  of  the  induced 
circulation  in  the  tropics. 
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Figure  8.  The  meridional  divergence  of  the  moisture  (10'  2  g/kg/s)  averaged  over  days  30-365, 
obtained  from  the  NSEAM  aquaplanet  experiments  with  (a)  L20e,  (b)  L30e,  and  (c)  L30u  vertical  grid. 


[22]  On  the  other  hand,  from  the  Hovmoeller  plot  it  is 
found  that  the  Kelvin  wave  propagation  speed  for  L30u 
(Figure  9c)  is  higher  than  that  for  L30e  (Figure  9b)  and  is 
closer  to  the  observed  values  while  it  is  similar  to  that  for 
L20e  (Figure  9a).  The  difference  in  the  vertical  scale 
inferred  from  the  cloud  top  distribution  (Figure  10)  can  be 
linked  again  to  this  speed  difference.  The  greater/smaller 
vertical  scale  for  L30u/L30e  correspond  to  the  higher/lower 
Kelvin  wave  propagation  speed  for  L30u  (Figures  9c)/L30e 
(Figure  9b). 


[23]  As  shown  in  this  study,  altering  the  vertical  structure 
of  the  model  systematically  modifies  the  moisture  budget  of 
the  boundary  layer,  characteristics  of  the  convective  precip¬ 
itation,  and  distribution  of  the  heating  in  clouds,  thereby 
influencing  the  equivalent  depth  of  the  modeled  atmosphere 
and  propagation  characteristics  of  the  equatorial  waves.  As 
in  the  case  with  respect  to  L20e  (Figure  7a),  the  smaller 
vertical  grid  intervals  of  L30e  are  associated  with  less 
moisture  than  L30u  (Figure  11a),  consistently  with  Inness 
et  al.  [2001],  Despite  the  less  moisture  present  in  the 
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Figure  9.  Time-longitude  plot  from  days  30  to  120  of  the  convective  precipitation  (10-1  mm/d) 
averaged  for  an  equatorial  band  (5°S-5°N)  obtained  from  the  NSEAM  aquaplanet  experiments  (a)  L20e, 
(b)  L30e,  (c)  L30u,  (d)  L30u  with  new  LCL,  and  (e)  L30u  with  modified  ATmax.  Note  that  a  smaller 
contour  interval  is  used  for  Figures  9c-9e.  The  solid  lines  are  drawn  to  help  estimate  the  propagation 
speed  of  cloud  clusters.  The  dashed  lines  in  Figures  9d  and  9e  are  drawn  to  identify  slower  moving 
clusters  distinguished  from  faster  moving  counterparts  denoted  by  the  solid  lines. 


boundary  layer,  however,  the  smaller  grid  intervals  of  L30e 
facilitate  easier  saturation  of  the  layers  and  thus  stronger 
precipitation  than  L30u  (Figure  6)  as  in  comparison  with 
L20e.  In  the  two  comparisons  (between  L30e  and  L20e  and 
between  L30e  and  L30u),  the  more  abundant  moisture 
content  in  L20e  and  L30u  may  merely  be  a  consequence 
of  less  precipitation  than  L30e. 

3.4.  Sensitivity  to  Model  Precipitation  Physics  and 
Simulation  of  Kelvin  Waves  and  Madden-Julian 
Oscillation 

[24]  The  spectral  analysis  of  tropical  waves  introduced  by 
Wheeler  and  Kiladis  [1999]  is  an  effective  and  convenient 
tool  for  investigating  the  organization  and  propagation  of 
equatorial  waves.  Figure  12  compares  the  (equator)  sym¬ 
metric  component  of  the  wave  number-frequency  decom¬ 
position  divided  by  the  background  spectrum  obtained 
using  the  simulated  convective  precipitation.  The  straight 
lines  in  Figure  12  denote  the  theoretical  shallow  water 
dispersion  relations  for  linear  Kelvin  waves  obtained  with¬ 
out  background  flow  and  with  the  equivalent  depths  of  50  m, 
25  m  and  12  m  based  on  observations  and  corresponding 


phase  speeds  of  22  m/s,  16  m/s  and  11  m/s  (Figure  12a; 
respectively  from  left).  The  relatively  low  propagation 
speeds  of  L20e  and  L30e  shown  by  the  Flovmoeller  plots 
of  the  convective  precipitation  (Figures  9a  and  9b)  are 
confirmed  in  Figures  12a  and  12b. 

[25]  Here,  we  investigate  to  what  extent  the  Kelvin  wave 
propagation  depends  on  details  of  the  model  physics  that 
involve  vertical  motion.  First,  we  modify  the  Emanuel 
convective  parameterization  so  that  that  a  twice-thicker 
layer  is  used  to  determine  the  LCL.  The  Hovmoeller  plot 
from  this  experiment  with  L30u  is  shown  in  Figure  9d.  The 
rather  scattered  cloud  clusters  that  were  present  in  L30u 
with  the  old  LCL  formulation  (Figure  9c)  are  now  more 
organized  with  the  new  LCL  formulation  (Figure  9d)  while 
maintaining  the  dominant  propagation  speed  of  the  old  LCL 
case  (Figure  9c). 

[26]  Next,  we  modify  the  value  of  the  specified  temper¬ 
ature  deficit  at  the  LCL  (ATmax)  in  the  Emanuel  convective 
cloud  parameterization  from  0.9  to  1.2  K.  This  value 
accounts  for  the  subgrid  turbulence  within  the  boundary 
layer  and  its  increase  would  enhance  the  chance  of  convec¬ 
tion.  The  result  (Figure  9e)  displays  far  stronger  organiza- 
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Figure  9.  (continued) 


tion  of  the  cloud  clusters  with  this  change.  It  is  interesting 
that  increasing  ATmax  amplifies  the  bimodal  structure  of  the 
eastward  propagating  perturbations  with  the  slow  moving 
waves  propagating  approximately  with  the  speed  of  MJO 
(marked  by  dashed  lines  in  Figure  9e)  and  faster  waves 
corresponding  to  observed  Kelvin  waves  (solid  lines). 
When  the  precipitation  fields  are  filtered  for  the  observed 
MJO  and  Kelvin  wave  modes  (Figure  13)  following 
Wheeler  and  Kiladis  [1999],  we  can  see  the  resemblance 
between  the  observed  Kelvin  and  MJO  modes  (Figure  13a) 
and  the  simulated  modified  ATmax  case  (Figure  13b),  with 
more  intense  Kelvin  waves  developing  during  the  “wet” 


MJO  phase  (red  contours  in  Figure  13).  This  relationship  is 
not  evident  in  other  experiments  presented  in  this  study. 

[27]  The  Wheeler-Kiladis  spectral  diagrams  of  the  con¬ 
vective  precipitation  corresponding  to  these  modified  phys¬ 
ics  experiments  are  presented  in  Figures  12d  and  12e.  In 
Figures  12d  and  12e,  two  modes  can  be  observed:  one  faster 
mode  corresponding  to  Kelvin  waves  with  period  of  about 
5-10  days  and  phase  speeds  of  about  1 0  m/s  and  the  other 
slower  mode  with  low  zonal  wave  numbers  and  periods 
between  30  and  50  days  corresponding  to  observed  MJO 
propagation  characteristics.  This  is  clearly  seen  as  separate 
spectral  peaks  particularly  in  Figure  12e  and,  additionally, 
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Figure  10.  Vertical  distribution  of  convective  cloud  top 
percentage  for  various  aquaplanet  simulations  (see  the  text 
for  details)  for  an  equatorial  band  (5°S-5°N)  averaged  over 
days  30-365.  The  convective  cloud  top  in  terms  of  the 
pressure  in  each  model  grid  box  was  counted  every  12  h. 
The  ordinate  (hPa)  is  in  logarithmic  scale,  roughly 
corresponding  to  the  elevation. 


in  Figure  12f  that  is  from  another  experiment  performed 
with  both  the  modified  LCL  and  ATmax. 

[28]  The  spectral  peaks  with  periods  greater  than  30  days 
at  eastward  propagating  wave  numbers  1  and  2  are  generally 
regarded  as  a  signature  of  the  MJO  [Neale  and  Hoskins, 
2001b]  while  higher  wave  numbers  (3-5)  are  sometimes 
implicitly  considered  to  fall  in  the  MJO  range  (e.g.,  Wheeler 
and  Kiladis  [1999,  Figure  3b]  based  on  estimated  outgoing 
longwave  radiation,  Masunaga  et  al.  [2006,  Figure  le] 
based  on  observed  deep  storm  cloudiness,  and  Cho  et  al. 
[2004,  Figure  5d]  based  on  observed  rainfall).  Our  spectrum 
(Figure  12f,  in  particular),  although  based  on  the  idealized 
aquaplanet  conditions,  is  fairly  similar  to  these  results.  In 
this  regard,  the  slower  moving  mode,  which  is  represented 
by  the  slower  moving  clusters  (Figures  9d  and  9e)  and  the 
low-frequency  spectral  peaks  (Figures  12c —  12f)  may  be 
associated  with  the  MJO. 

[29]  We  attempt  to  further  investigate  the  strong  spectral 
peak  that  appears  at  wave  number  5  in  the  L30u_LCL2  case 
(Figure  12d),  which  is  also  given  by  Neale  and  Hoskins 
[2001b,  Figure  3b]  and  Masunaga  et  al.  [2006,  Figure  le]. 
From  Figure  10  we  find  that  the  L30u_LCL2  case  contains 
two  cloud  tops  of  similar  percentage  (200  hPa  and  300  hPa) 
whereas  L30u  contains  the  major  top  at  300  hPa.  The  higher 
and  lower  cloud  tops  of  L30u_LCL2  may  be  associated 
with  the  faster  (near  the  theoretical  lines)  and  slower 
(around  50  day  period  and  5E  wave  number)  wave  modes, 
respectively  (Figure  1 2d).  The  higher  cloud  tops  of 
L30u_LCL2  is  taller  than  that  ofL30u,  but  the  corresponding 
vertical  scale  is  not  necessarily  deeper  because  of  the 
secondary  (lower)  top  at  300  hPa,  which  is  at  the  same 
height  as  the  major  top  of  L30u  and  thus  makes  the  two 
vertical  scales  comparable.  This  may  be  why  the  two  cases 
produce  relatively  similar  overall  spectral  distributions  as 
depicted  in  Figures  12c  and  12d.  It  is  also  noted,  although 
not  very  evident,  that  the  L30u_LCL2  case  has  slightly 


higher  percentage  for  the  lower  top  at  300  hPa  than  the 
L30u_LCL2  +  dTm2  case  whereas  it  has  lower  percentage 
for  the  higher  top  at  200  hPa.  The  more  dominant  (i.e., 
higher  percentage)  lower  top  at  300  hPa  for  L30u_LCL2 
may  then  lead  to  a  shallower  depth  scale,  than  that  would  be 
otherwise  without  the  lower  top,  which  may  have  reduced 
the  wave  propagation  speed  (not  shown).  In  this  regard,  the 
strong  spectral  peak  of  L30u_LCL2  (Figure  12d)  at  5E  and 
around  50  day  period  may  be  a  result  of  a  shallower  depth 
scale  which  pushed  the  spectrum  further  toward  the  lower- 
frequency  range  and  may  thus  represent  a  very  slow  Kelvin 
wave  that  is  in  the  MJO  frequency  range. 

[30]  However,  we  would  like  to  note  here  that  our 
analysis  based  on  the  vertical  depth  scale  of  the  convective 
processes,  which  is  inferred  from  the  vertical  cloud  top 
pressure  percentage  distribution,  is  limited  and  may  not 
provide  the  full  picture,  especially  because  the  cloud  top 
pressures  do  not  reflect  the  existence  of  shallower  clouds 
under  the  deepest  cloud  within  a  given  model  column.  As 


Average  over  700-950  hPa 


Figure  11.  (a)  The  difference  between  L30e  and  L30u  of 
the  water  vapor  (10“'  g/kg)  and  (b)  the  comparison 
between  L30e  and  L30u  of  the  cumulus-induced  heating  rate 
(10  1  K/d)  derived  from  the  temperature  tendency  by  the 
Emanuel  convective  cloud  parameterization.  The  data  have 
been  averaged  between  700  and  950  hPa  over  days  30-365 
of  the  NSEAM  aquaplanet  simulations. 
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Figure  12.  The  symmetric  component  of  the  wave  number-frequency  decomposition  of  the  convective 
precipitation  following  Wheeler  and  Kiladis  [1999]  for  cases  (a)  L20e,  (b)  L30e,  (c)  L30u,  (d)  L30u  with 
new  LCL,  (e)  L30u  with  modified  ATmax,  and  (f)  L30u  with  both  new  LCL  and  modified  ATmax  for  days 
30  through  365.  See  the  text  for  the  descriptions  of  the  straight  lines  and  the  speed  values. 


Masunaga  et  al.  [2006]  pointed  out,  the  high-wave  number 
MJO  mode  may  simply  have  been  exaggerated  because  of 
the  normalization  by  the  arbitrarily  defined  background 
spectrum  and  contribute  to  highly  dispersive  nature  of  the 
MJO  [see  also  Sal  by  and  Hendon,  1994], 

4.  Conclusions  and  Further  Remarks 

[31]  Even  under  a  controlled  environment  of  the  aqua- 
planet  experiments,  the  propagation  of  equatorial  atmo¬ 


spheric  waves  vary  widely  among  models;  even  the 
simulated  propagation  directions  of  the  equatorial  waves 
are  opposite  in  some  models  [see,  e.g.,  Moncrieff  et  al., 
2007,  Figure  3].  This  can  be  due  to  many  reasons,  but  the 
convective  cloud  parameterization  is  the  prime  suspect. 
Atmospheric  models  are  sensitive  to  many  details  of  the 
model,  in  particular  to  the  physical  parameterizations  in¬ 
volving  vertical  motions  of  air  parcels.  Indeed,  cumulus 
parameterization  is  considered  one  of  the  most  dominant 
modeling  components  that  generate  uncertainties  in  large- 
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(a)  Kelvin  waves  Sc  MJO  (obs) 
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(b)  Kelvin  waves  Sc  MJO  (NSEAM) 
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Figure  13.  The  powers  in  the  convective  precipitation  ((mm/h)2)  of  Kelvin  wave  amplitude  (color 
shades)  and  wet  (i.e.,  active)  phase  of  MJO  (red  contours),  obtained  by  filtering  the  precipitation  data  in 
time  and  space  [Wheeler  and  Kiladis,  1999]:  (a)  2006  TRMM  (Tropical  Rainfall  Measuring  Mission) 
observation  and  (b)  NSEAM  L30u  aquaplanet  simulation  with  modified  ATmax.  In  typical  NSEAM 
simulations,  MJO  and  Kelvin  waves  are  not  necessarily  correlated  unlike  the  observation  shown  in 
Figure  13  a,  but  the  observed  relationship  starts  appearing  when  ATmax  is  increased  to  make  it  easier  to 
trigger  convection  as  shown  in  Figure  13b.  Note  the  color  shade  values  for  Figure  13b  are  multiplied  by 
10  to  facilitate  the  comparison  with  those  for  Figure  13a;  that  is,  the  color  bar  values  for  Figure  13b  range 
from  0.8  to  4. 


scale  atmospheric  simulation  [e.g.,  Tomita  et  al,  2005; 
Moncrieff  et  al.,  2007]  (see  also  Lee  et  al.  [2003]  who 
investigated  the  sensitivity  of  aquaplanet  simulations  to 
different  cumulus  parameterizations,  which  implicitly  dem¬ 
onstrates  the  uncertainties). 

[32]  In  this  study  we  have  explored  such  uncertainties 
utilizing  our  new  prototype  global  atmospheric  forecast 
model,  NSEAM.  We  have  further  demonstrated  that  the 
simulations  of  convective  activity  can  be  quite  different 
even  with  the  same  convective  cloud  parameterization 
depending  on  the  distribution  of  vertical  levels.  We  have 
also  demonstrated  that  the  simulations  are  sensitive  to  the 
details  of  the  parameterization  itself  such  as  the  LCL-related 
processes,  which  modify  the  organization  and  propagation 
of  equatorial  waves.  For  this  reason,  we  suggest  some 
“standardization”  of  vertical  levels  for  more  direct  inter¬ 
comparison  among  the  models  (e.g.,  for  the  APE).  While 
we  have  presented  results  suggesting  that  the  new  model 
can  be  calibrated  to  capture  some  basic  characteristics  of  the 
Kelvin  waves  and  MJO,  the  present  study  only  provides 
limited  analyses  of  the  reasons  for  the  sensitivity  of  the 
simulations  to  various  model  parameters,  and  mainly  points 


out  the  seriousness  of  the  uncertainties  concerning  convec¬ 
tive  processes  and  their  coupling  to  large-scale  wave  motion. 

[33]  Although  not  reported  in  this  study,  in  order  to 
further  investigate  the  effect  of  the  variation  of  the  vertical 
resolution,  we  performed  additional  experiments  by  sepa¬ 
rately  adding  1  or  2  more  levels  in  L30e  and  removing 
1  level  from  L30u  at  low  levels  and  also  by  separately 
adding  1  or  2  more  levels  in  L30u  around  <7—0.85.  The 
results  were  qualitatively  similar  to  the  original  cases  (i.e., 
L30e  and  L30u),  implying  that  our  aquaplanet  simulation  is 
not  highly  sensitive  to  a  local  addition  or  removal  of  one  or 
two  levels.  A  nonlocal,  systematic  modification  of  the 
layers,  however,  induces  systematic  response  of  the  model 
as  has  been  shown  in  this  study.  We  investigated  the 
dependence  of  the  aquaplanet  simulations  on  the  integration 
time  step  ( Williamson  and  Olson  [2003]  performed  a 
comprehensive  sensitivity  study  on  this  subject).  In  addition 
to  the  default  time  step  of  300  s,  we  also  tried  150  s  and 
450  s  for  L20e  cases.  We  found  from  the  corresponding 
Wheeler-Kiladis  diagrams  that  the  phase  speed  was  propor¬ 
tional  to  the  time  step  to  some  extent  with  the  150  s  case 
yielding  the  fastest  and  closest  to  the  theoretical  shallow 
water  values.  We  also  performed  a  lower  horizontal  resolu- 
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tion  counterpart  (3  elements  and  8  polynomials)  of  the 
aquaplanet  experiments.  Despite  the  lower  resolution,  the 
simulation  is  fairly  similar  to  the  higher-resolution  counter¬ 
part,  except,  for  example,  that  the  precipitation  band  over  the 
equator  is  wider  and  less  homogeneous.  While  Lorant  and 
Royer  [2001]  reported  more  intensified  and  concentrated 
convective  cells  in  their  higher-resolution  GCM  simulations, 
our  preliminary  results  show  less  difference  due  possibly  to 
the  strong  controlled  environment  of  the  aquaplanet  condi¬ 
tions  and  the  selected  magnitude  of  the  viscosity. 

[34]  In  conclusion,  the  results  presented  in  this  study 
demonstrate  that  the  tropical  convective  precipitation  phys¬ 
ics  and  associated  equatorial  wave  dynamics  are  highly 
sensitive  to  details  of  the  vertical  structure  as  well  as 
parameterized  physical  processes  of  the  model.  The  observed 
sensitivity  may  explain  the  difficulties  in  modeling  of  the 
equatorial  waves  and  MJO  in  large-scale  models  of 
the  atmosphere.  There  exists  complex  interaction  between 
the  dynamics  and  physics  in  terms  of  the  circulation  around 
the  convective  area  along  the  equator,  which  is  likely  to  be 
underrepresented  in  the  models.  Further  investigation  is 
needed  to  obtain  more  fundamental  and  insightful  under¬ 
standing  of  the  sensitivity  and  the  interaction  found  in  this 
study. 


Appendix  A:  Derivation  of  the  Second-Order 
Horizontal  Viscosity  Operator  and  Its  Extension  to 
High-Order  Hyperviscosity 

[35]  A  K  =  1  hyperviscosity  operator  given  in  equation 
(4)  can  be  constructed  in  the  following  way  using  integra¬ 
tion  by  parts,  based  on  Giraldo  [1999]; 

J  1 pV2q  dQ  —  j  ip nV q  dT  —  j  VipVqdQ  (Al) 
n  r  n 


where  ip  represents  the  basis  function,  n  the  normal  vector, 
Q  refers  to  the  interior  domain  and  T  to  the  boundaries.  For 
a  global  model,  periodicity  is  the  only  lateral  boundary 
condition  and,  thereby,  forces  the  first  term  on  the  right- 
hand  side  to  vanish.  Therefore,  at  discrete  grid  points  “1” 
we  can  write  the  Laplacian  operator  as  follows 

V2?i  =  -Mr'  Lijqj  (A2) 

where  “j”  denotes  neighboring  points  and 

Mi  =  Wi\Ji\,  Lij  =  Wi\Ji\\7ipij  Vipjj  (A3) 

are  the  mass  and  Laplacian  matrices,  respectively,  and  J  is 
the  Jacobian  with  corresponding  quadrature  weight,  w.  One 
can  now  extend  this  idea  further  to  construct  an  arbitrarily 
high-order  hyperviscosity  operator  as  follows 


V  qt  = 


n  i~Mi{  l‘j)> 


?/• 


(A4) 
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